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TWO SPECTRAL METHODS FOR 2D QUASI-PERIODIC 
SCATTERING PROBLEMS 


KUI DU* 

Abstract. We consider the 2D quasi-periodic scattering problem in optics, which has been 
modelled by a boundary value problem governed by Helmholtz equation with transparent boundary 
conditions. A spectral collocation method and a tensor product spectral method are proposed to 
numerically solve the problem on rectangles. The discretization parameters can be adaptively chosen 
so that the numerical solution approximates the exact solution to a high accuracy. Our methods also 
apply to solve general partial differential equations in two space dimensions, one of which is periodic. 
Numerical examples are presented to illustrate the accuracy and efficiency of our methods. 
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1. Introduction. With the development of electromagnetic and optics technol¬ 
ogy and the increasing demands in industrial and military applications, scattering 
from periodic structures |29| have attracted much interest in recent years. A variety 
of numerical methods including finite difference methods m, finite element meth¬ 
ods [l[5l[T0l[Ti[71|9l|40l[22], spectral and spectral element methods [271 US] : integral 
equation methods [iiKniESKn], Dirichlet-to-Neumann map methods |28j , and mode 
expansion method |3Ql[T] have been developed by the engineering community and the 
applied mathematical community for solving linear diffraction problems from periodic 
structures. 

Finite difference and finite element methods are easy to implement and result in 
sparse linear systems, which enable the use of sparse direct solvers. However, these 
methods typically have significant dispersion errors for high-frequency problems and 
thus fail to provide accurate solutions. If the medium is piecewise constant, boundary 
integral equations formulated on the interfaces of multilayer structures m or the 
boundaries of the obstacles m M are natural and mathematically rigorous. By 
exploiting high-order quadratures, one obtains much higher efficiency and accuracy 
than finite difference and finite element methods. For problems with general medium, 
spectral methods [2311381 da daiM] can be employed. They are simple to implement 
and typically require relatively small number of unknowns to attain a fixed accuracy. 
We refer to |35j for a brief survey of existing spectral methods for partial differential 
equations. 

In this paper, we consider 2D quasi-periodic scattering problems. By introducing 
transparent boundary conditions, the problem is governed by a Helmholtz equation 
with variable coefficients defined on rectangular domains. We propose a spectral col¬ 
location method and a tensor product spectral method. The first one uses the spectral 
collocation techniques [23l [38] and the second one, based on separable representations 
of the differential operator, combines the Fourier spectral method and the ultraspher- 
ical spectral method [32] . For vertically layered medium case, both methods only 
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require to solve a one dimensional problem. The two spectral methods proposed in 
this paper can adaptively determine the number of unknowns and approximate the 
solution to a high accuracy. We remark that the tensor product spectral method leads 
to matrices that are banded or almost banded. By employing the fast algorithm in 
[32) . a layered medium approximation problem can be used as a preconditioner for 
the problem with general media. 

Recently, Chebfun software [20] was extended to solve problems in two dimen¬ 
sions [361135] . both of which are nonperiodic. We refer to [39] for the summary of 
the mathematics and algorithms of Chebyshev technology for nonperiodic functions. 
Extension to periodic functions in one dimension was proposed in [42| . Our spectral 
methods can be used to solve the problems in two space dimensions, one of which is 
periodic. Implementing our methods with the adaptive strategies of Chebfun software 
is easy. 

The rest of this paper is organized as follows. In §2 the model problem is for¬ 
mulated and further reduced to a boundary value problem. In §3 we present the 
spectral collocation method for the problem. Section 4 is devoted to the tensor prod¬ 
uct spectral method. Numerical examples illustrating the accuracy and efficiency of 
the methods are reported in §5. We present brief concluding remarks in §6. 

2. Problem setting. Assume that no currents are present and that the fields 
are source free. Then the electromagnetic fields in the whole space are governed by 
the following time-harmonic Maxwell equations (time dependence e““^‘) 

(2.1) V X E = iw/rH, 

(2.2) V X H = -iweE, 

where i = \/—1 is the imaginary unit, uj is the angular frequency, jj, is the permeability, 
£ is the permittivity, E is the electric field and H is the magnetic field. 

In this paper, we consider a simple quasi-periodic scattering model: /r is constant 
everywhere; e is periodic in x variable of period 27r, invariant in z variable, and 
constant away from the region |y| < 1 — d, i.e., there exist constants £“*■ and £“ such 
that 


£{x,y,z) = e+, for ?/ > 1 - 5, 

£{x,y,z) = £~, for y< 5—1. 

In the TM polarization, the electric field E takes the simpler form 

E(x, y, z) = E(a;, y) = (0, 0, u(x, y)). 

The Maxwell equations (I2.1D - (I2.2D yield the Helmholtz equation: 

(2.3) Au + io^£yu = 0. 

Consider the plane wave u' = jg incident from the above, where ao = 

w\/e+/rsin0, /3o = uj^/£+ y, cos 9 , and —7r/2 < 9 < 7r/2 is the angle of incidence with 
respect to the positive y-axis. The incident wave it' leads to reflected wave it'' and 
transmitted wave it*. For y > 1, we have it = it' -I- itand for y < — 1, it = it*. We are 
interested in quasi-periodic solution u with phase e'“‘’^'', i.e., ite“'““''' is 27r-periodic 
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in X variable. Therefore, the reflected and transmitted waves can be written as 

= , y>l, 

jGZ 

= y<-l, 

jGZ 

where Vj and tj are unknown complex scalar coefficients and 
aj = ao + j, 

/3j = yJuj'^e+fj.- a'^, Im(/3j) > 0, 

7j = Im(7j) > 0. 

Note that if u)^J£^y, is real, then the j’s satisfying fij > 0 correspond propagating 
modes. Throughout, we assume that (3j ^ 0 and 7 j ^ 0 for all j £ Z. This assumption 
excludes the “resonant” cases, where waves can propagate along the i-axis. 

For a quasi-periodic function /(a;), define the linear operators S and T by [8l[3l|6] 

(2.4) (,5/)(x)=^i^,/,e-^^ 

jez 

(2.5) (r/)(a;) = ^i7,/,e'“^", 

jez 

where 

a'ez 

1 

/, = ^ I /(x)e-‘“^Mx. 

Let V denote the unit outward normal. We obtain transparent boundary conditions: 

(2.6) = —2i/3oe'““^“'^'’, on (0,27r)x{l}, 

(2.7) d^u-Tu = ^, on (0,27r) x {-!}. 

The quasi-periodic scattering problem is to solve the Helmholtz equation (12.31) in 
the rectangular domain H := (0, 27r) x (—1,1) subject to the transparent boundary 
conditions (I2.6I) - ()2.7I) and the quasi-periodic boundary condition 

(2.8) it(x + 27r, 2 /) = e‘“°^’^'u(a:, 2 /). 

Define v = ite“‘“°“^. Then v satisfies 


AagV + ui’^eyv = 0, 

in 

L! = (0,2^)x (-1,1), 

dyv - = -2i/3oe-'^«, 

on 

(0,2^) x{l}, 

dyv + = 0, 

on 

(0,2^)x{-l}, 

v{x + 2TT,y) =v{x,y), 

in 



where the operator is defined by = dxx + — |q;oP + dyy The solution 

to (12.91) is unique at all but a discrete set of frequencies uj when the incident angle 9 
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is fixed. Existence and uniqueness of the solution is strictly proved in Dobson [TS] by 
a variational approach. 

Since the operators S and T given in (I2.4I) - (I2.5I) are defined by infinite series, the 
computation has to be truncated in practice. For simplicity, we truncate S and T as 
follows: 


N-q 

N-q 

{Nf){x)= Y 

where is a sufficiently large integer and q is an integer satisfying that \ < q < N 
and all propagating modes are contained in the middle of the truncated series. Next, 
we propose two spectral methods for the following problem 


( 2 . 10 ) 


AtjoU -I- uj'^eyv = 0, 

in 

n = i0,27r) X (-1,1), 

dyv - = -2i/3oe-'^°, 

on 

(0,2^) x{l}. 

dyv + e-i“o^r^(e“»^u) = 0, 

on 

(0,2^)x{-l}, 

vlx + 2TT,y) = v{x,y), 

in 



We refer to [40l Section 3.1] for the discussion on the existence and uniqueness of the 
solution of the problem (12.101) . In this paper, we focus on the accurate and efficient 
spectral methods for the problem (12.101) . Thus we always assume that the discrete 
problem has a unique solution. 

3. A spectral collocation method. We approximates the problem (12.101) by 
Fourier discretization in x variable and Chebyshev discretization in y variable. Let 
Xn = nh and ym = cos{rmr/M) with h = 2tt/N , n = 1,2,..., N, and m = 0,1,..., M. 
Introduce the first-order Fourier differentiation matrix [33|1T1[3H], 


N odd , 

N even , 

the second-order Fourier 


(F = 


(F = 


0, 

* = 


* 3, 

o' 

i=3, 


i^3, 


differentiation matrix [25l|4ll[38], D,, = 


N odd , 


N even , 


1 

r ^2 

1 




1 3^2 

+ 12’ , 

“ 3)d 
-^cot 
2 

(*-J> 

2 

s 

[ (-1)*- 

-^■+1 ^ CSC ■ 

2 



1 



— < 


“ 6’ 

-3)h 

2 ’ 

1=3, 

a^j ^ 

[ (-1)*- 

-^■+1 ^ csc^ 

2 



* = j, 

* ^ 3, 
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and the (M+1) x (M+1) Chebyshev differentiation matrix [I11|1T1|3H], T)y = 


II 

o 

2M^ + 1 

- 

^MM — 

2M2 + 1 

6 

6 

dl = 

Vi 

* = 1,. 

..,M-1, 

2(1-2/f)’ 

II 

Q (-1)*+^ 

cj {y^ -VoV 


bJ =0,1 

Ci = < 

f 2, i = 0 or M, 

[ 1, otherwise. 



The discretization of the Helmholtz equation in (12.101) takes the form 

(3.1) PVDj^ + 2iaoPVDT - |aopPV + PD^V + uJ‘^^lP{e © V) = 0, 

where P is the (M — 1) x (M + 1) matrix obtained by deleting the hrst and last rows of 
the identity matrix Im+i, V = [vmn]mLo n=iJ '^mn is the approximate solution at the 
point {Xn, Vm), £ = [£mn]m=o,n=i^ ^rnn = £{xn, Um), and 0 denotes the componentwise 
product. The operators e“'““^S'^(e*““^-) and e“*“o®T^(e*“““*) can be approximated 
by TV X iV matrices: 

(3.2) S = G*A; 3 G := G*diag{i^i_„ i^ 2 -„ • • • , i^A-JG, 

(3.3) T = G*A^G := G*diag{i7i_„ i72-„ • • • , i7A-JG, 

where G is the N x N matrix with {i,j) entry /^/N, i,j = 

The discrete forms of the transparent boundary conditions are given by 

(3.4) ejDj^V - eJVS'^ = -2i^oe“'^“e'^, 

(3.5) elByV + VT^ = 0, 

where 

eo = [1 0 ••• 0]^ G R^+\ 
eM = [0 ••• 0 1]^ G 
e = [1 1 ••• l]"^ G K^. 

Let 

(3.6) H = (D^^ + 2iaoDx - |q!oPIa) © P + Ia © (PDy) + w^^(Ia © P)diag{vec(e)}. 
Combining (13.11) . (13.41) and (13.51) yields the linear system 

(3.7) Av = g, 
where 

Ia © (ejDy) - S © ej 
A = H 

Ia © (sAflly) + T © 
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and 


V = vec(V), 


g = 


—2i/3oe ‘^°e 

0 

0 


g f^MN+N 


Remark 3.1. The matrix P in (EU is called the downsampling matrix [m 
obtained by using the inner Chebyshev points of the second kind, 



as resampling points. One might choose other points, such as the Chebyshev points of 
the first kind, 

/ (2m — l)7r\ 

Vm = cos - 1) j ’ m = 1,2,...,M-1. 

We refer to |43j for explieit construction of rectangular differentiation matrix PDy. 

3.1. Layered medium case. Assume that the medium in 12 is vertically lay¬ 
ered, i.e., e{x,y) = e{y). Then, the matrix H in (j3.6|) takes the form 

H = {'Da:x + 2iaoDa; - |aoplAr) 0 P -|- Iat 0 (P(D^ -|- 

where 


= diag{e(?/o), e(yi), • • • , e(2/M)}. 

Let F be the discrete Fourier transform matrix with {i,j) entry /n j^ 

i, j = 1,2,... ,N. We have the following propositions. The proofs are straightforward. 
Proposition 3.2. The matrices S in (|3^ and T in dSEl) satisfy 

S = F*AsF, T = F*AtF, 

where 

As = diag{Af,A^,--- ,A^} 

= diag{i^o,i^i,--- , i^i-g, i^2-g, • • • 

and 


At = diag{A7,A],--- ,Xjf} 

= diag{i7o,i7i,--- , i7Ar_q, i7i_q, i72-q, • • • ,i7-i}- 

Proposition 3.3. The first-order spectral differentiation matrix D^, and the 
second-order spectral differentiation matrix Gxx can be diagonalized by F. Specifically, 
we have 


D, = F*A,F, D,, = F*A,,F, 

where 

A, =diag{A?,A^,-- - ,A^} 

^ r i • diag{0,1, • • • , (A - l)/2, (1 - N)/2, (3 - N)/2, ■ ■ ■ , -1}, if N is odd, 
\ i • diag{0,1, • • • , A/2 -1,0,1- A/2, 2 - A/2, • • • , -1}, if N is even. 
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A,, = diag{Ar,Ar,--- ,A^"} 

^ f -diag{0,1,4, • • • , ((TV - l)/2)2, ((1 - N)/2)'^, • • • , 4,1}, if N is odd, 
\ -diag{0,1,4, • • • , {N/2f, {N/2 - 1)^, • • • , 4,1}, if N is even. 


Multiplying the linear system (13.711 by 


from the left yields 
(3.8) 
where 

V = (F 0 Im+i)v 


F 0 Im-1 


In ® (ejDy) - As 0 ej 
+ 2 iaoA 2 ; - laoplA) 0 P + Ijv 0 (P(D2 + 
In 0 {eljBy) + At 0 


V = g, 


g = 


—2i/?oe 

0 

0 


e = Fe 


y/N 

0 


Reordering the unknowns and equations in (13.81) . we obtain N uncoupled linear sys¬ 
tems of order M -|- 1: 


ejDy - i/3oeJ 


■ -2i/3oe-'^“y]V ■ 

P(D2+a;V^)- laopP 

Vl = 

0 

e^Dy -P iyoe^ 


0 


and for z = 2,..., A, 


ejD - Af ej 

(Af" + 2iaoAf - |aoP)P V P(D2 + uj^fiey) 

eJ^Dy + X]eJ^ 


Vi = 0, 


where Vi is the z-th column of VF"*". Obviously, Vi = 0 for z 1. After solving the 
linear system (13.9L we obtain v by v = (F* 0 Im-i-i)v. 

4. A tensor product spectral method. We propose a tensor product spectral 
method (see §4.3) for the problem (12.101) by combining the Fourier spectral method 
and the ultraspherical spectral method [35] • 

4.1. Fourier spectral method for periodic second order linear ordinary 
differential equations. We consider the second-order linear ordinary differential 
equation (ODE) 


(4.1) Cw := a{x)w"{x) -I- b{x)w'{x) + c{x)w{x) = f{x), x C [0, 27r] 


with periodic boundary conditions. Here, a{x), b{x), c{x) and f{x) are periodic 
functions on [0, 27r]. The Fourier spectral method finds an infinite vector 

I'r 

W-2 W-i Wq Wi W2 ■ ■ ■ , 


w = 
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such that the Fourier expansion of the solution of (HU) is given by 

w(x) = X G [0,27r]. 


Note that we have 


w'(x) = w"(x) = y^(— 

jez jez 


Wje 


.ya: 


The first-order differentiation operator is given by 


= i 


-3 


-2 


-1 


and the second-order differentiation operator is given by "D^. In order to handle vari¬ 
able coefficients of the forms a{x)w"{x), b{x)w'{x) and c{x)w{x) in (14.11) . we need to 
represent the multiplication of two Fourier series as an operator on coefficients. Let 
T[a] denote the multiplication operator that represents multiplication of two Fourier 
series, i.e., if w is a vector of Fourier expansion coefficients of w(x), then T[a]w re¬ 
turns the Fourier expansion coefficients of a{x)w{x). Suppose that a{x) is given by 
its Fourier series 


a{x) = 

jez 

Then the explicit formula for T[a] is given by the following Toeplitz operator 



r[a] = 


a-2 

a_i 

tto 

ai 

0-2 


This multiplication operator looks dense; however, if a(x) is approximated by a 
trigonometric polynomial of degree n, then T[a] is banded with a bandwidth of n. 
Combining the differentiation and multiplication operators yields 

(r[ap2 + -f r[c])w = f, 
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where w and f are vectors of Fourier expansion coefficients of w(x) and f(x), respec¬ 
tively. We need truncate the operator to derive a practical numerical scheme. Let Pjv 
be the projection operator satisfying 

PjVW = [ 0 Iat 0 ]w = [ Wl-g W2-q ■■■ WN-q ] "^ ■ 

We obtain the following linear system 

{VnC°VI){Vn^) = 


where 


= r[a]vi + r[b]v, + r[c]. 

The truncation parameter N can be adaptively chosen so that the numerical solution 
approximates the exact solution to relative machine precision. 

4.2. The ultraspherical spectral method |32| for second order linear 
ordinary differential equations. Consider the second order linear ordinary differ¬ 
ential equation 

(4.2) Cw := a{x)w''{x) + b{x)w'{x) + c{x)w{x) = /(x), 

where a(x), b{x), c(x), w{x) and f{x) are functions defined on [—1,1]. The ultras¬ 
pherical spectral method finds an infinite vector 

w = [ Wo ixi • ■ ■ ] 

such that the Chebyshev expansion of the solution of (14.21) is given by 

OO 

w{x) ='^WjTj{x), xG[-1,1], 

i=o 

where Tj (x) is the degree j Chebyshev polynomial. 

Note that we have the following recurrence relation 

[ 2"-in(A-l)!Cii\, n > A, 

dx^ lO, O^n^A — 1, 


where is the ultraspherical polynomial with an integer parameter A > 1 of degree 
J EH- Then the differentiation operator for the Ath derivative is given by 


(4.3) 


Vx = 2^-i(A- 1)! 


0 A 

A+1 

A-f 2 


A > 1. 


Here 0 in (|4.3I) denotes a A-dimensional zero row vector. The conversion operator 
converting a vector of Chebyshev expansion coefficients to a vector of expansion 
coefficients, denoted by Sq, and the conversion operator converting a vector of 
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expansion coefficients to a vector of expansion coefficients, denoted by Sx, are 

given by 


■ 1 0 

1 

2 




■ 1 0 

A 

A+2 



1 

0 

1 



A 

0 

_^ 


2 

2 



A-l-1 

A+3 



1 

2 

0 


II 


A 

A+2 

0 




1 





A 




2 





A-l-3 


- 










We also require the multiplication operator Alo[o] that represents multiplication of 
two Chebyshev series, and the multiplication operator AlA[a] that represents multi¬ 
plication of two series, i.e., if w is a vector of Chebyshev expansion coefficients 
of w{x), then Alo[a]w returns the Chebyshev expansion coefficients of a{x)w{x), and 
Af ■ ■ ■ 5ow returns the expansion coefficients of a[x)w[x). Suppose that 

a{x) has the Chebyshev expansion 


a{x) = 

3=0 


Then Alo[a] can be written as [35]: 


2 ao 

ai 

02 

03 ••• ■ 


■ 0 

0 

0 

0 • 


ai 

2ao 

Ol 

02 

1 

+ 2 

Ol 

02 

03 

04 


02 

Ol 

2oo 

ai 

02 

03 

Q4 

O5 


03 

02 

Ol 

2ao 


03 

04 

O5 

06 











- 


The explicit formula for the entries of AIaIo] with A > 1 is given in [35]. This mul¬ 
tiplication operators M.x{<A with A > 0 look dense; however, if a(x) is approximated 
by a truncation of its Chebyshev or series, then AIaM is banded. Combining 
the differentiation, conversion and multiplication operators yields 

(4.4) [A42[o.]'D2 + iSiAf-4- iSiiSoAlo[c])w = iSiiSof, 

where w and f are vectors of Chebyshev expansion coefficients of w{x) and f{x), 
respectively. We need truncate the operator to derive a practical numerical scheme. 
Let Qm be the projection operator given by 

Qm = [ Im 0 ] . 

We obtain the following linear system 

(Qm/:®Qm)(Qmw) = {QMSiSoQli){QMi), 

where 

£» = M2[a]V2 + SiMi[b]Vi + SiSoMo[c]. 
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4.3. Discretization of (|2.10|) for the special case e{x,y) = 4>{x)ip{y). We 
seek the approximate solution v{x,y) to the problem (12.101) . 

M-l N-q 

i=0 j=l-q 

The linear differential operator in (12.101) can be expressed as 

+ 1® + uj'^yLcf) ® 

where £“ := dxx + 2iao9x — |aop! Cy := X is the identity operator, and (p and ip 
are the corresponding multiplication operations. The discretization of the Helmholtz 
equation in (12.101) takes the form 

(4.5) QCVX^ + QYV + w = 0, 

where Q is the {M — 2) x M matrix obtained by deleting the last two rows of the 
identity matrix 1m, C = Qm^hSqQ^, X = T’Ni'Dl + 2iaQVx)Vjj - |aoplAr, Y = 
QMV 2 Qlt, ^ = VnUcPYPI, ^ = Qm5i5oMoMQ^, V = The 

discrete forms of the transparent boundary conditions are given by 

(4.6) b^V - a^VA^ = -2i/?oe-‘^°eT 

(4.7) bJV + aJVA.y = 0, 

where 

Kp = diag{i/3i_g, i^ 2 - 9 , • • • , 

A^ = diag{i 7 i_,, i72-g, • • • , i7A-g}, 

bi = [0 1 4 ••• {M-lff, 

b2 = [0 1 -4 ••• (-l)^(M- 1)2]T, 

ai = [1 1 ••• 1]^ G 

a2 = [1 -1 ••• (-1)(^-^)]^ G 


and Gq is the g-th column of the identity matrix I^r. 

Let 

(4.8) H = X 0 (QC) + Iat 0 (QY) + ® (Q’®'). 

Combining (14.51) . (14.61) and (14.71) yields the linear system 


(4.9) 

where 


and 


Av = g, 


A = 


Iat 0 bj" - A/3 0 aj" 

Iat 0 bj + A.y 0 aj , 

H 


g = 


-2i/3oe“ 

0 

0 


G C 


MN 


V = vec(V) 
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Remark 4.1. For general bivariate function e{x,y), we can use its low rank 
approximant, i.e., sum of functions of the form where (f^x) and ip{y) are 

univariate functions. An algorithm [Ml 133 which is mathematically equivalent to 
Gaussian elimination with complete pivoting can he used to construct low rank ap¬ 
proximations. 


4.3.1. Layered medium case. In this case, we have e{x, y) = i.e., 4>{x) = 

1. Then = I^r. The matrix H in (14.81) takes the form 

H = X (g) (QC) + Iat g) (QY + wVQ’®')- 


Reordering the unknowns and equations in (14.9L we obtain 



b?" - i/^oa?" 


■ -2i/3oe-*'5'> ■ 

(4.10) 

bj -f iqoaj 

V, = 

0 


_ Q(Y + a;2y'®'-|aopC) _ 


0 


and 

(4.11) 


bj + hj-q^2 

Q(Y - {{j - qf + 2ao{j - q) + laoHC) 


Vj = 0, 


j 7^ 9, 


where v^- is the j-th column of V. Obviously, Vj = 0 for j ^ q. 

Note that the fast algorithm in |M] can be used to solve the linear system (14.101) . 
which requires operations. Here, m is the number of Chebyshev points 

needed to resolve the function ijj{y). The coefficient matrix resulting from a layered 
medium approximation problem can be used as a preconditioner for the problem with 
general medium. The corresponding computational complexity for the preconditioner 
solve is 0{m?MN). 


5. Numerical results. We have performed numerical experiments in numerous 
cases. In this section, we present a few typical results of these experiments to illus¬ 
trate the accuracy and efficiency of the two spectral methods. All computations are 
performed with MATLAB R2012a. 

The parameters are chosen as 0 = Stt/T, o; = 10, £“*■ = = I, ^ = I. We 

consider three cases: 


r+4 


ei{x,y) = e{y) = I -f e^^’ 

^ ^^g^+4-ycos(.sinW2))^ 


The medium characterized by ei(x, y) is vertically layered. The bivariate function 
£ 2 ( 2 :, y) is of rank 2. We use the rank 7 approximant obtained by the algorithm 
proposed in 131132 ] to approximate e^{x,y). In Figured we plot the three media 
and the real parts of the corresponding waves obtained by our spectral methods. We 
observe that the numerical results obtained by the two methods coincide well for all 
the three media. 

We also test the performance of the problem with £i(a:,y) as a preconditioner 
for the problem with e 3 (x,y). The fast algorithm in [32] is used to solve the linear 
systems (14.101) and (14.111) . The GMRES algorithm [33] is used as the iterative solver. 
The initial guess is set to be the zero vector. GMRES with the layered medium 
preconditioner obtained a solution with the relative residual norm less than 10“® 
at the 37th iteration, while GMRES without preconditioning almost stagnates; see 
Figured for the convergence history. 
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Fig. 1. Media (Top) and corresponding real parts of the waves obtained by the spectral col¬ 
location method (Middle) and by the tensor product spectral method (Bottom). From left to right: 
numerical results for £i{x,y), £ 2 (x,y), and £ 3 {x,y), respectively. 



Fig. 2. Convergence history for GMRES applied to the linear system 114.911 . 

6. Concluding remarks. We have proposed a spectral collocation method and 
a tensor product spectral method for solving the 2D quasi-periodic scattering prob¬ 
lem. Both of the methods can adaptively determine the number of unknowns and 
approximate the solution to a high accuracy. The tensor product spectral method 
is more interesting because it leads to matrices that are banded or almost banded, 
which enable the use of the fast, stable direct solver [52] ■ Based on this fast solver, a 
layered medium preconditioning technique is used to solve the problem with general 
media. Our methods also apply to solve general partial differential equations in two 
space dimensions, one of which is periodic. Extension of our methods to bi-periodic 
structure diffraction grating problem |3] is being considered. 


REFERENCES 


[1] B. Bai and J. Turunen, Fourier modal method for the analysis of second-harmonic generation 
in two-dimensionally periodic structures containing anisotropic materials, J. Opt. Soc. 












14 


KUI DU 


Amer. B, 24 (2007), pp. 1105-1112. 

[2] G. Bao, Finite element approximation of time harmonic waves in periodic structures, SIAM 

J. Numer. Anal., 32 (1995), pp. 1155-1169. 

[3] -, Numerical analysis of diffraction by periodic structures: TM polarization, Numer. 

Math., 75 (1996), pp. 1-16. 

[4] -, Variational approximation of Maxwell’s equations in biperiodic structures, SIAM J. 

Appl. Math., 57 (1997), pp. 364-381. 

[5] G. Bao, Y. Gao, and H. Yang, Numerical solution of diffraction problems by a least-squares 

finite element method. Math. Methods Appl. Sci., 23 (2000), pp. 1073-1092. 

[6] G. Bao and Y. Chen, A nonlinear grating problem in diffractive optics, SIAM J. Math. Anal., 

28 (1997), pp. 322-337. 

[7] G. Bao, Z. Chen, and H. Wu, Adaptive finite-element method for diffraction gratings.. Journal 

of the Optical Society of America. A, Optics, image science, and vision, 22 (2005), pp. 1106— 
1114. 

[8] G. Bao, D. 0. Dobson, and J. A. Cox, Mathematical studies in rigorous grating theory, J. 

Opt. Soc. Amer. A, 12 (1995), pp. 1029-1042. 

[9] G. Bao, Y. Li, and H. Wu, Numerical solution of nonlinear diffraction problems, J. Comput. 

Appl. Math., 190 (2006), pp. 170-189. 

[10] G. Bao and H. Yang, A least-squares finite element analysis for diffraction problems, SIAM 

J. Numer. Anal., 37 (2000), pp. 665—682 (electronic). 

[11] A. Barnett and L. Greengard, A new integral representation for quasi-periodic fields and its 

application to two-dimensional band structure calculations, J. Comput. Phys., 229 (2010), 
pp. 6898-6914. 

[12] - , A new integral representation for quasi-periodic scattering problems in two dimensions, 

BIT, 51 (2011), pp. 67-90. 

[13] J. P. Boyd, Chebyshev and Fourier spectral methods, Dover Publications, Inc., Mineola, NY, 

second ed., 2001. 

[14] 0. Canuto, M. Y. Hussaini, a. Quarteroni, and T. a. Zang, Spectral methods in fluid 

dynamics. Springer Series in Computational Physics, Springer-Verlag, New York, 1988. 

[15] C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang, Spectral methods: Funda¬ 

mentals in single domains. Scientific Computation, Springer-Verlag, Berlin, 2006. 

[16] Z. Chen and H. Wu, An adaptive finite element method with perfectly matched absorbing 

layers for the wave scattering by periodic structures, SIAM J. Numer. Anal., 41 (2003), 
pp. 799-826. 

[17] M. H. Cho and a. H. Barnett, Robust fast direct integral equation solver for quasi-periodic 

scattering problems with a large number of layers. Optics express, 23 (2015), pp. 1775-1799. 

[18] D. C. Dobson, Optimal design of periodic antireffective structures for the helmholtz equation, 

European Journal of Applied Mathematics, 4 (1993), pp. 321-339. 

[19] T. A. Drisgoll and N. Hale, Rectangular spectral collocation, IMA Journal of Numerical 

Analysis, to appear (2015). 

[20] T. A. Driscoll, N. Hale, and L. N. Trefethen, Chebfun guide, Pafnuty Publications, Ox¬ 

ford, 2014. 

[21] K. Du AND W. Sun, A finite difference method for second harmonic generation from periodic 

structures, manuscript, (2015). 

[22] K. Du AND M. Zhang, A tensor product finite element method for the diffraction grating 

problem with transparent boundary conditions, submitted, (2015). 

[23] B. Fornberg, a practical guide to pseudo spectral methods, vol. 1 of Cambridge Monographs on 

Applied and Computational Mathematics, Cambridge University Press, Cambridge, 1996. 

[24] A. CiLLMAN AND A. Barnett, A fast direct solver for quasi-periodic scattering problems, J. 

Comput. Phys., 248 (2013), pp. 309-322. 

[25] D. Gottlieb, M. Y. Hussaini, and S. A. Orszag, Theory and applications of spectral meth¬ 

ods, in Spectral methods for partial differential equations (Hampton, Va., 1982), SIAM, 
Philadelphia, PA, 1984, pp. 1-54. 

[26] Y. He, M. Min, and D. P. Nicholls, A spectral element method with transparent boundary 

condition for periodic layered media scattering, to appear, (2015). 

[27] Y. He, D. P. Nicholls, and J. Shen, An efficient and stable spectral method for electromag¬ 

netic scattering from a layered periodic structure, J. Comput. Phys., 231 (2012), pp. 3007- 
3022. 

[28] Y. Huang and Y. Y. Lu, Scattering from periodic arrays of cylinders by Dirichlet-to-Neumann 

maps. Journal of lightwave technology, 24 (2006), pp. 3448-3453. 

[29] J. D. JOANNOPOULOS, S. C. JOHNSON, J. N. WiNN, AND R. D. Meade, Photonic Crystals: 

Molding the Flow of Light, Princeton University Press, second ed., 2008. 





TWO SPECTRAL METHODS FOR QUASI-PERIODIC SCATTERING 


15 


[30] B. Maes, P. Bienstman, and R. Baets, Modeling second-harmonic generation by use of mode 

expansion, J. Opt. Soc. Amer. B, 22 (2005), pp. 1378—1383. 

[31] F. W. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, NIST handbook of mathe¬ 

matical functions, Cambridge University Press, 2010. 

[32] S. Olver and A. Townsend, A fast and well-conditioned spectral method, SIAM Rev., 55 

(2013), pp. 462-489. 

[33] Y. Saad and M. H. Schultz, GMRES: a generalized minimal residual algorithm for solving 

nonsymmetric linear systems, SIAM J. Sci. Statist. Comput., 7 (1986), pp. 856-869. 

[34] J. Shen, T. Tang, and L.-L. Wang, Spectral methods, vol. 41 of Springer Series in Computa¬ 

tional Mathematics, Springer, Heidelberg, 2011. Algorithms, analysis and applications. 

[35] A. Townsend and S. Olver, The automatic solution of partial differential equations using a 

global spectral method, J. Comput. Phys., to appear (2015). 

[36] A. Townsend and L. N. Trefethen, An extension of Chebfun to two dimensions, SIAM J. 

Sci. Comput., 35 (2013), pp. C495—C518. 

[37] A. Townsend and L. N. Trefethen, Gaussian elimination as an iterative algorithm, SIAM 

News, 46 (2013). 

[38] L. N. Trefethen, Spectral methods in MATLAB, vol. 10 of Software, Environments, and Tools, 

Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2000. 

[39] -, Approximation theory and approximation practice. Society for Industrial and Applied 

Mathematics (SIAM), Philadelphia, PA, 2013. 

[40] Z. Wang, G. Bao, J. Li, P. Li, and H. Wu, An adaptive finite element method for the 

diffraction grating problem with transparent boundary conditions, SIAM J. Numer. Anal., 
53 (2015), pp. 1585-1607. 

[41] J. A. C. Weideman and S. C. Reddy, A MATLAB differentiation matrix suite, ACM Trans. 

Math. Software, 26 (2000), pp. 465—519. 

[42] G. B. Wright, M. Javed, H. Montanelli, and L. N. Trefethen, Extension of chebfun to 

periodic functions, SIAM J. Sci. Comput. submitted, (2015). 

[43] K. Xu AND N. Hale, Explicit construction of rectangular differentiation matrices, IMA Journal 

of Numerical Analysis, to appear (2015). 



